#ifndef _GHOSTFILLER_H_
#define _GHOSTFILLER_H_

#include "RegularGrid/RectDomain.h"
#include "RegularGrid/Func.h"
#include "Core/Tensor.h"


template <int Dim> class GhostFiller;

template <int Dim>
class GhostFiller{
public:
  using iVec = Vec<int,2>;
  using rVec = Vec<Real,2>;
  
  GhostFiller(const RectDomain<Dim>& aDomain);
  
  template <class TFunc>
  void fillAllSides(Tensor<Real,Dim>& res, const char* BCType, const
                    TFunc& func);

  void fillAllSides(Tensor<Real,Dim>& res, const char* BCType, int);
protected:
  template <class TFunc>
  void fillOneSide(Tensor<Real,Dim>& res, int dim, int side, const char
                   BCType, const TFunc& func);

  void fillOneSide(Tensor<Real,Dim>& res, int dim, int side, const char
                   BCType, int);
protected:
  RectDomain<Dim> Domain;
};

template <int Dim>
GhostFiller<Dim>::GhostFiller(const RectDomain<Dim>& aDomain){
  assert(aDomain.getCentering() == NodeCentered);
  assert(aDomain.getNumGhost() == 1);
  Domain = aDomain;
}

// Case Dim = 2 -----------------------------------------------

template<>
template <class TFunc>
void GhostFiller<2>::fillOneSide(Tensor<Real,2>& res, int dim, int
                                 side, const char BCType, const TFunc& func){
  const Box<2>& bx = Domain;
  const iVec& lo = bx.lo();
  const iVec& sz = bx.size();
  const rVec& dx = Domain.spacing(); 
  if (dim == 1 && side == -1){
     if (BCType == 'D'){
       for(int i = 0 ; i < sz[0] ; i++){
         iVec Node{i,0};
         rVec rNode = (Node - lo)*dx;
         res(0,i) = func.F(rNode);
         res(-1,i) = 6*res(1,i)-8*res(2,i)+3*res(3,i);
       }
     }
     else if (BCType == 'X'){
       for(int i = 0 ; i < sz[0] ; i++){
         res(-1,i) = 6*res(1,i)-8*res(2,i)+3*res(3,i);
       }
     }
     else if (BCType == 'N'){
       for(int i = 0 ; i < sz[0] ; i++){
         iVec Node{i,0};
         rVec rNode = (Node - lo)*dx;
         res(-1,i) = res(1,i)-2*dx[1]*func.Fy(rNode);
       }
     }
  }
  else if (dim == 1 && side == 1){
     if (BCType == 'D'){
       for(int i = 0 ; i < sz[0] ; i++){
         iVec Node{i,sz[1]-1};
         rVec rNode = (Node - lo)*dx;
         res(sz[1]-1,i) = func.F(rNode);
         res(sz[1],i) = 6*res(sz[1]-2,i)-8*res(sz[1]-3,i)+3*res(sz[1]-4,i);
       }
     }
     else if (BCType == 'X'){
       for(int i = 0 ; i < sz[0] ; i++){
         res(sz[1],i) = 6*res(sz[1]-2,i)-8*res(sz[1]-3,i)+3*res(sz[1]-4,i);
       }
     }
     else if (BCType == 'N'){
       for(int i = 0 ; i < sz[0] ; i++){
         iVec Node{i,sz[1]-1};
         rVec rNode = (Node - lo)*dx;
         res(sz[1],i) = res(sz[1]-2,i)+2*dx[1]*func.Fy(rNode);
       }
     }
  }
  else if (dim == 0 && side == -1){
     if (BCType == 'D'){
       for(int i = 0 ; i < sz[1] ; i++){
         iVec Node{0,i};
         rVec rNode = (Node - lo)*dx;
         res(i,0) = func.F(rNode);
         res(i,-1) = 6*res(i,1)-8*res(i,2)+3*res(i,3);
       }
     }
     else if (BCType == 'X'){
       for(int i = 0 ; i < sz[1] ; i++){
         res(i,-1) = 6*res(i,1)-8*res(i,2)+3*res(i,3);
       }
     }
     else if (BCType == 'N'){
       for(int i = 0 ; i < sz[1] ; i++){
         iVec Node{0,i};
         rVec rNode = (Node - lo)*dx;
         res(i,-1) = res(i,1)-2*dx[0]*func.Fx(rNode);
       }
     }
  }
  else if (dim == 0 && side == 1){
     if (BCType == 'D'){
       for(int i = 0 ; i < sz[1] ; i++){
         iVec Node{sz[0]-1,i};
         rVec rNode = (Node - lo)*dx;
         res(i,sz[0]-1) = func.F(rNode);
         res(i,sz[0]) = 6*res(i,sz[0]-2)-8*res(i,sz[0]-3)+3*res(i,sz[0]-4);
       }
     }
     else if (BCType == 'X'){
       for(int i = 0 ; i < sz[1] ; i++){
         res(i,sz[0]) = 6*res(i,sz[0]-2)-8*res(i,sz[0]-3)+3*res(i,sz[0]-4);
       }
     }
     else if (BCType == 'N'){
       for(int i = 0 ; i < sz[1] ; i++){
         iVec Node{sz[0]-1,i};
         rVec rNode = (Node - lo)*dx;
         res(i,sz[0]) = res(i,sz[0]-2)+2*dx[0]*func.Fx(rNode);
       }
     }
  }
  else
    assert(0);
}

template<>
void GhostFiller<2>::fillOneSide(Tensor<Real,2>& res, int dim, int side, const char
                                 BCType, int){
  D2Func0 f;
  fillOneSide(res,dim,side,BCType,f);
}

template<>
template <class TFunc>
void GhostFiller<2>::fillAllSides(Tensor<Real,2>& res, const char*
                                  BCType, const TFunc& func){
  for(int i = 0 ; i < 4 ; i++){
    if (BCType[i] == 'N')
      fillOneSide(res,i/2,-1+(i%2)*2,BCType[i],func);
  }
  for(int i = 0 ; i < 4 ; i++){
    if (BCType[i] != 'N')
      fillOneSide(res,i/2,-1+(i%2)*2,BCType[i],func);
  }
}

template<>
void GhostFiller<2>::fillAllSides(Tensor<Real,2>& res, const char*
                                  BCType, int){
  for(int i = 0 ; i < 4 ; i++){
    if (BCType[i] == 'N')
      fillOneSide(res,i/2,-1+(i%2)*2,BCType[i],0);
  }
  for(int i = 0 ; i < 4 ; i++){
    if (BCType[i] != 'N')
      fillOneSide(res,i/2,-1+(i%2)*2,BCType[i],0);
  }
}

// Case Dim = 1 -----------------------------------------------

template<>
template<class TFunc>
void GhostFiller<1>::fillOneSide(Tensor<Real,1>& res, int dim, int
                                 side, const char BCType, const
                                 TFunc& func){
  assert(dim == 0);
  const Box<1>& bx = Domain;
  const Vec<int,1>& lo = bx.lo();
  const Vec<int,1>& sz = bx.size();
  const Vec<Real,1>& dx = Domain.spacing(); 
  if (side == -1){
    if (BCType == 'D'){
      Vec<int,1> Node{0};
      Vec<Real,1> rNode = (Node - lo)*dx;
      res(0) = func.F(rNode);
      res(-1) = 6*res(1)-8*res(2)+3*res(3);
    }
    else if (BCType == 'X'){
      res(-1) = 6*res(1)-8*res(2)+3*res(3);
    }
    else if (BCType == 'N'){
      for(int i = 0 ; i < sz[0] ; i++){
        Vec<int,1> Node{0};
        Vec<Real,1> rNode = (Node - lo)*dx;
        res(-1) = res(1)-2*dx[0]*func.Fx(rNode);
      }
    }
  }
  else if (side == 1){
    if (BCType == 'D'){
      Vec<int,1> Node{sz[0]-1};
      Vec<Real,1> rNode = (Node - lo)*dx;
      res(sz[0]-1) = func.F(rNode);
      res(sz[0]) = 6*res(sz[0]-2)-8*res(sz[0]-3)+3*res(sz[0]-4);
    }
    else if (BCType == 'X'){
      res(sz[0]) = 6*res(sz[0]-2)-8*res(sz[0]-3)+3*res(sz[0]-4);
    }
    else if (BCType == 'N'){
      Vec<int,1> Node{sz[0]-1};
      Vec<Real,1> rNode = (Node - lo)*dx;
      res(sz[0]) = res(sz[0]-2)+2*dx[0]*func.Fx(rNode);
     }
  }
  else
    assert(0);
}

template<>
void GhostFiller<1>::fillOneSide(Tensor<Real,1>& res, int dim, int side, const char
                                 BCType, int){
  D1Func0 f;
  fillOneSide(res,dim,side,BCType,f);
}

template<>
template<class TFunc>
void GhostFiller<1>::fillAllSides(Tensor<Real,1>& res, const char*
                                  BCType, const TFunc& func){
   fillOneSide(res,0,-1,BCType[0],func);
   fillOneSide(res,0,1,BCType[1],func);
}

template<>
void GhostFiller<1>::fillAllSides(Tensor<Real,1>& res, const char*
                                  BCType, int){
  fillOneSide(res,0,-1,BCType[0],0);
  fillOneSide(res,0,1,BCType[1],0);
}
#endif //_GHOSTFILLER_H_
